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Astract. In this paper, we study the Local Discontinuous Galerkin methods for 
nonlinear, time-dependent convection- diffusion systems. These methods are an ex- 
tension of the Runge-Kutta Discontinuous Galerkin methods for purely hyperbolic 
systems to convection-diffusion systems and share with those methods their high 
parallelizability, their high-order formal accuracy, and their easy handling of com- 
plicated geometries, for convection dominated problems. It is proven that for scalar 
equations, the Local Discontinuous Galerkin methods are L 2 -stable in the nonlinear 
case. Moreover, in the linear case, it is shown that if polynomials of degree k are 
used, the methods are fc-th order accurate for general triangulations; although this 
order of convergence is suboptimal, it is sharp for the LDG methods. Preliminary 
numerical examples displaying the performance of the method are shown. 
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1. Introduction. In this paper, we study the Local Discontinuous Galerkin 
(LDG) methods for nonlinear, convection-diffusion systems of the form 

3,u + V-F(u,Du) = 0, in (0, T) x ft, 

where ft C and u = (mj, . . . . u m ) # . The LDG methods are ail extension of the 
Runge-Ivutta Discontinuous Galerkin (RIvDG) methods for the nonlinear hyper- 
bolic system 

dfVL + V * f(u) = 0, in (0. T) x ft, 

introduced by the authors [12,13,14,15,16], and further developed by Atkins and 
Shu [2], Bassi and Rebay [4], Bey and Oden [7], Biswas. Devine, and Flaherty [8], 
deCougny et al [17], Devine et al [19, 20], Lowrie, Roe and van Leer [30], and by 
Oztiiran et al [33]. The RIvDG methods are constructed by applying the explicit 
time discretizations introduced by Shu [37] and Shu and Osher [38,39] to a space 
discretization that uses discontinuous basis functions. Since the space discretization 
is highly local in character and produces easily invertible, block-diagonal mass 
matrices and since the time-marching scheme is explicit, the RIvDG method is a 
highly parallelizable method; see Biswas, Devine, and Flaherty [8]. Moreover, it is 
not only a formally high-order accurate method that can easily handle complicated 
geometries, but it satisfies a cell entropy inequality that enforces a nonlinear L 2 - 
stability property even without the slope limiters typical of this method; see Jiang 
and Shu [27]. 

Extensions of the RIvDG method to hydrodynamic models for semiconductor 
device simulation have been constructed by Chen et al [9], [10]. In these extensions, 
approximations of the derivatives of the discontinuous approximate solution are 
obtained by using a simple projection into suitable finite elements spaces. This 
projection requires the inversion of global mass matrices, which in [9] and [10] are 
Tumped’ in order to maintain the high parallelizabilitv of the method. Since in [9] 
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and [10] polynomials of degree one are used, the ‘mass lumping’ is justified; however, 
if polynomials of higher degree were used, the ‘mass lumping’ needed to enforce 
the full parallelizability of the method could cause a degradation of the formal 
order of accuracy. Fortunately, this is not an issue with the methods proposed by 
Bassi and Rebay [5] (see also Bassi et al [6]) for the compressible Navier-Stokes 
equations. In these methods, the original idea of the RIvDG method is applied 
to both u and D u which are now considered as independent unknowns. Like the 
RIvDG methods, the resulting methods are highly parallelizable methods of liigh- 
order accuracy which are very efficient for time-dependent, convection-dominated 
flows. The LDG methods are a generalization of these methods. 

The basic idea to construct the LDG methods is to suitably rewrite the convection- 
diffusion system into a larger, degenerate, first-order system and then discretize it 
by the RIvDG method. By a careful choice of this rewriting, nonlinear stability can 
be achieved even without slope limiters, just as the RIvDG method in the purely 
hyperbolic case; see Jiang and Shu [27], In the linear case, the stability result leads 
to the sub-optimal rate of (Aj)‘ : for the L°°(0, T;L 2 )-norm of the error if polynomi- 
als of degree at most k are used. However, these estimates are sharp, as numerical 
evidence reported in Bassi et al. [6] and in this paper indicate. In the purely 
hyperbolic case, the rate of convergence of ( Ar) fc+1 / 2 is recovered, as expected. 
Indeed, this is the same rate of convergence obtained for the original Discontinuous 
Galerkin method (introduced by Reed and Hill [35]) for purely hyperbolic case by 
Johnson and Pitkaranta [28] and confirmed to be optimal by Peterson [34]. LeSaint 
and Raviart [29] proved a rate of convergence of (Ax) k for general triangulations 
and of (Aj')k +1 for Cartesian grids; Richter [36] obtained the optimal rate of con- 
vergence of ( Aa-)^ +1 for some structured two-dimensional non-Cartesian grids. The 
technique for proving the error estimates used in this paper is different from the 
techniques used in the above mentioned papers. It is very simple and relies, as ex- 
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pected, on a straightforward combination of (i) the Instability of the LDG method 
and of (ii) the approximation properties of the finite element spaces. 

The LDG methods introduced in this paper are very different from the so-called 
Discontinuous Galerkin (DG) method for parabolic problems introduced by Jamet 
[26] and studied by Eriksson, Johnson, and Thomee [25], Eriksson and Johnson [21, 
22, 23, 24], and more recently by Makridakis and Babuska [31]. In the DG method, 
the approximate solution is discontinuous only in time, not in space; in fact, the 
space discretization is the standard Galerkin discretization with continuous finite 
elements. This is in strong contrast with the space discretizations of the LDG 
methods which use discontinuous finite elements. To emphasize this difference, we 
call the methods developed in this paper the Local Discontinuous Galerkin meth- 
ods. We also must emphasize that the large number of degrees of freedom and the 
restrictive conditions of the size of the time step for explicit time-discretizations, 
render the LDG methods inefficient for diffusion- dominated problems; in this sit- 
uation, the use of methods with continuous-in-space approximate solutions is rec- 
ommended. However, as for the successful RKDG methods for purely hyperbolic 
problems, the extremely local domain of dependency of the LDG methods allows 
a very efficient parallelization that by far compensates for the extra number of 
degrees of freedom in the case of convection-dominated flows. 

Many researchers have worked in the devising and analysis of numerical meth- 
ods for convection-dominated flows. In particular, Dawson [18] and, more recently, 
Arhogast and Wheeler [1] have developed and analyzed methods that share several 
properties with the LDG methods: They use discontinuous-in-space approxima- 
tions, they are locally conservative, and they approximate the diffusive fluxes with 
independent variables (by using a mixed method). We refer the reader interested in 
numerical methods for convection-dominated flows to [18] and [1] and the references 
therein. 
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Another numerical method that uses discontinuous approximations is the one 
proposed and studied by Baker et al. [3]. This method, however, is not for 
convection- dominated flows but for the Stokes problem. The advantage of using 
discontinuous approximations in this setting is that this allows for a pointwise ver- 
ification of the incompressibility condition at the interior of each triangle. Optimal 
estimates are obtained. 

In this paper, we restrict ourselves to the semidiscrete LDG methods for convection- 
diffusion problems with periodic boundary conditions. Our aim is to clearly display 
the most distinctive features of the LDG methods in a setting as simple as possible. 
The fully discrete methods for convection-diffusion problems in bounded domains 
will be treated in a forthcoming paper. This paper is organized as follows: In §2, 
we introduce the LDG methods for the simple one-dimensional case d — 1 in which 

F(u,Du) — f(u) — a(n)d x u, 

u is a scalar and a(u) > 0 and show some preliminary numerical results displaying 
the performance of the method. In this simple setting, the main ideas of how to 
devise the method and how to analyze it can be clearly displayed in a simple way. 
Thus, the L 2 -stability of the method is proven in the general nonlinear case and 
the rate of convergence of (Ax)* in the L°°(0, T;L 2 )-norm for polynomials of degree 
k > 0 in the linear case is obtained; this estimate is sharp. In §3, we extend these 
results to the case in which u is a scalar and 

F, ( u , Du) = f< O ) - a t _,(«) d JCj u , 

\< } <d 

where a,y defines a positive semidefinite matrix. Again, the L 2 -stability of the 
method is proven for the general nonlinear case and the rate of convergence of 
(Ax)* in the L°°(0, T;L 2 )-norm for polynomials of degree k > 0 and arbitrary tri- 
angulations is proven in the linear case. In this case, the multidimensionality of the 
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problem and the arbitrariness of the grids increase the technicality of the analysis 
of the method which, nevertheless, uses the same ideas of the one-dimensional case. 
In §4, the extension of the LDG method to multidimensional systems is briefly 
described and concluding remarks are offered. 

2. The LDG methods for the one-dimensional case. In this section, we 
present and analyze the LDG methods for the following simple model problem: 

dt u + d j: (f(u) — a(u)d r v) = 0 in (0, T) x (0, 1), (2.1a) 

u(t = 0) = i/o, on (0, 1). (2.1b) 

with periodic boundary conditions. 

a. General formulation and main properties. To define the LDG method, 
we introduce the new variable q = \Jo( u ) d r u and rewrite the problem (2.1) as 
follows: 


di u + d x (/(«) - yja(u)q) = 0 

i—l 

O 

X 

h 

0 

• rH 

(2.2a) 

q - Or g(u) = 0 

in (0. T) x (0. 1), 

(2.2b) 

u{t = 0 ) = 1/0. 

on (0, 1). 

(2.2c) 


where g(u) = f u \J a(s) ds . The LDG method for (2.1) is now obtained by simply 
discretizing the above system with the Discontinuous Galerkin method. 

To do that, we follow [13] and [14]. We define the flux h = ( h u Ji q ) / as follows: 

h (u,q) = ( /(«) - y/a(u)q, -g(u) )'. (2.3) 

For each partition of the interval (0, 1). { fj+ 1/2 we I j = i x j-i/ 2 - x j+i/ 2 )' 

and A.ij = x ;+1 / 2 — .r | /■> for j = 1. .... A' ; we denote the quantity maxi^xw Ax j 
l)v Ax . We seek an approximation w* = (uh,qhY to w = ( u ,q)' such that for 
each time / £ [0, T], both </*(#) and qi, {t) belong to the finite dimensional space 

v h = V h k = {v € I'(0, 1) : v\ h € j = 1 N}. 


(2.4) 
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where P k {I) denotes the space of polynomials in I of degree at most k. In order 
to determine the approximate solution («/,,<//,), we first note that by multiplying 
(2.2a), (2.2b), and (2.2c) by arbitrary, smooth functions v u , v q , and u t , respectively, 
and integrating over Ij, we get, after a simple formal integration by parts in (2.2a) 
and (2.2b), 


/ dt u{xj)v u (x)dx - f h u (w(xj))d x v u (x) dx 
Jlj Jij 

+ h u (w(x j+1/2 ,t)) v u (x~ +1/2 ) - h u (w{x J _ l/2 J))v u (x+_ 1/2 ) = 0. 
I q(x,t) Vg(x) dx — f h q (w(xj))d x v q (x)dx 

Jlj Jlj 

+ h q (w{x j+l/2 J))v q (x~ +l/2 ) - h q (w{x } _ l/2 ,t))v q {x+_ 1/2 ) = 0. 
I u ( x , 0 ) Vi ( J" ) dx = I uq ( x ) Vi ( X ) dx . 

Jlj Jij 


(2.5a) 


(2.5b) 

(2.5c) 


Next, we replace the smooth functions v u , v q , and v, by test functions <>/, ,, and 
Vh,i, respectively, in the finite element space \\ and the exact solution w = {u,qY 
by the approximate solution w h = ( n /, - <lh ) 1 • Since this function is discontinuous 
in each of its components, we must also replace the nonlinear flux h(w(a 'j+i/ 2 , t )) 
by a numerical flux h(w) J+ 1 / 2 (f) = {h u {wf, ) j+l / 2 (t), h q {w h ) J+ 1 / 2 (t)) that will be 
suitably chosen later. Thus, the approximate solution given by the LDG method 
is defined as the solution of the following weak formulation: 


/ d t u h {x,t)v h ' U {x)dx - / h u {vr h (x,t))d x v htU (x)dx 

Jij Jij 

+ h u (w h ) j+1/2 (t)v htU {xj +l/2 ) - /»«(w A )>-i /2 (*) VA,.(xt_ 1/2 ) - 0, 

/ qh{x,t) Vh, q (x)dx - / h q (w h (x,t))d x v h , q (x)dx 
J i i Jii 

+ hq(*fh)j + l/2(t) Vh,q(*j +l/2 ) ~h q ( Wft )j-- [ / 2 {t)Vh, q {x+_ 1/2 ) = 0, 

J Uhi-r,0) Vh,i(x)dx — J u 0 (x) vh,i(x)dx, 


Vv h , u 6 P k (Ijh 

(2.6a) 


y v h , q e P k (Ij), 
(2.6b) 


Vt-i, e P k (ij). 


(2.6c) 



It only remains to choose the numerical flux h(w/, j /->(/). We use the notation: 

[p] = P + - P~< and p = 1( p + + p ~ ), 

and j /2 = p{-v ± +l /2 )• To be consistent with the type of numerical fluxes used in 
the RKDG methods, we consider numerical fluxes of the form 

h(w*) >+1/2 (#) = h(w fc (ar7 +1/2 ,i),w A (a-t + i /2 ,f)), 

that (i) are locally Lipschitz and consistent with the flux h, (ii) allow for a local 
resolution of qu in terms of u h , (iii) reduce to an E-flux (see Oslier [32]) when 
«(•) = 0, and that (iv) enforce the L 2 -stability of the method. 

To reflect the convection-diffusion nature of the problem under consideration, 
we write our numerical flux as the sum of a convective flux and a diffusive flux: 

h(w~, w + ) = h c „„„(w _ , w + ) + h rf , // (w _ .w + ). (2.7a) 

The convective flux is given by 

hrom-(w", W + ) = (/(«", « + ),0)', (2.7b) 

where /( </ “, m + ) is any locally Lipschitz E-flux consistent with the nonlinearity /, 
and the diffusive flux is given by 

hrf«//(w“, w+) = ( - T l' -g(v) y - c d t ff [w], (2.7c) 

where 


Cdiff = {-cn C o)' 

( 2 . 7 d ) 

C \2 = Ci 2 (w - , w+) is locally Lipschitz, 

(2.7e) 

0 
III 

C 

1 

in 

N 

c 

(2.7f) 


We claim that this flux satisfies the properties (i) to (iv). 
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Let us prove our claim. That the flux h is consistent with the flux h easily 
follows from their definitions, (2.3) and (2.7). That h is locally Lipschitz follows 
from the fact that /(•, •) is locally Lipschitz and from (2.7d); we assume that /(•) 
and a( ) are locally Lipschitz functions, of course. Property (i) is hence satisfied. 

That the approximate solution qh can be resolved element by element in terms of 
itf, by using (2.6b) follows from the fact that, by (2.7c), the flux h q - —g{u)—c \2 [ u ] 
is independent of qi, . Property (ii) is lienee satisfied. 

Property (iii) is also satisfied by (2.7f) and by the construction of the convective 
flux. 

To see that the property (iv) is satisfied, let us first rew r rite the flux h in the 
following way: 


f, _ +, /[</>(«)] [$(«)]_ — T n 

h(w ,w + ) = (- rn r — T — —</(«) ) -C[w], 


where 


c 


<Tl Cyz 

— Cl2 0 


cn 


i /[</>(«)] 


/(«.-, «+) 


MV M 

with <p{u) defined by 4>(u) = f" f(s)ds. Since /(•,•) is an E-flux, 


('ll = 


i r + 

\ / (/(-) 

m] Ju- 


f( u , ) ds > 0. 


and so, by (2.7d), the matrix C is semipositive definite. The property (iv) follows 
from this fact and from the following result. 


Proposition 2.1. (Instability) We have , 


\ [ n 2 h (x,T)dx + f f ql{xj)dxdf + 0 t , c ([ w />]) [ ul(x)dx 

- Jo Jo Jo - Jo 


where 


©r,d[w*]) = f JZ { [wa(#)]*C [w*(f )] | dt. 

l<j<N ^ ' j+1/2 
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This result will be proven in section §2.c. Thus, this shows that the flux h given 
by (2.7) does satisfy the properties (i) to (iv)- as claimed. 

Now, we turn to the question of the quality of the approximate solution defined 
by the LDG method. In the linear case f r = c and o(-) = a, from the above stability 
result and from the the approximation properties of the finite element space V/, , 
we can prove the following error estimate. We denote the L 2 (0. 1 )-norm of the 7-tli 
derivative of u by | u \(. 


Theorem 2.2. (L 2 -error estimate) Let e be the approximation error w — w/, . Then 
we have, 



| ( u (•*% T) |~ dx + 


ff 


( q (.rJ) | 2 dx dt + 0r,c( [ e ] 


1/2 

< C(Ax) k , 


where C — C(k , | a |*+ 1 , | u U+ 2 )' Ike purely hyperbolic case a — 0. the constant 
C is of order (A.r) 1 / 2 and in the purely parabolic case c — 0, the constant C is of 
order Ax for even values of k for uniform grids and for C identically zero . 


This result will be proven in section §2.d. The above error estimate gives a 
suboptimal order of convergence, but it is sharp for the LDG methods. Indeed. 
Bassi et al [6] report an order of convergence of order k + 1 for even values of 
k and of order k for odd values of k for a steady state, purely elliptic problem 
for uniform grids and for C identically zero. Our numerical results for a purely 
parabolic problem give the same conclusions; see Table 5 in the section §2.b. 

Our error estimate is also sharp in that the optimal order of convergence of 
k -fi 1/2 is recovered in the purely hyperbolic case, as expected. This improvement 
of the order of convergence is a reflection of the semipositive definiteness of the 
matrix C, which enhances the stability properties of the LDG method. Indeed, 
since in the purely hyperbolic case 
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the method enforces a control of the jumps of the variable uj, , as shown in Propo- 
sition 2.1. This additional control is reflected in the improvement of the order of 
accuracy from k in the general case to A’ + 1 /2 in the purely hyperbolic case. 

However, this can only happen in the purely hyperbolic case for the LDG meth- 
ods. Indeed, since c\\ = 0 for c — 0, the control of the jumps of Uh is not enforced 
in the purely parabolic case. As indicated by the numerical experiments of Bassi et 
al. [C] and those of section §2.b below, this can result in the effective degradation of 
the order of convergence. To remedy this situation, the control of the jumps of iif, 
in the purely parabolic case can be easily enforced by letting c 1 1 be strictly positive 
if | c | + | « | > 0. Unfortunately, this is not enough to guarantee an improvement 
of the accuracy: an additional control on the jumps of qh is required! This can be 
easily achieved by allowing the matrix C to be symmetric and positive definite when 
a > 0. In this case, the order of convergence of A- + 1/2 can be easily obtained for 
the general convection-diffusion case. However, this w r ould force the matrix entry 
('22 to be nonzero and the property (ii) of local resolvability of qh in terms of 
would not be satisfied anymore. As a consequence, the high parallelizability of the 
LDG would be lost. 

The above result shows how strongly the order of convergence of the LDG meth- 
ods depend on the choice of the matrix C. In fact, the numerical results of section 
§2.b in uniform grids indicate that with yet another choice of the matrix C. see 
(2.9), the LDG method converges with the optimal order of Ar + 1 in the general 
case. The analysis of this phenomenon constitutes the subject of ongoing work. 


b. Preliminary numerical results. 

In this section we provide preliminary numerical results for the schemes discussed 
in this paper. We will only provide results for the following one dimensional, linear 
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convection diffusion equation 


d t u + cd x u — a dj : u — 0 in (0, T) x (0, 2 zr), 
u(t = 0, x) = sin(,r ). on (0, 2 zr), 


where c and a > 0 are both constants; periodic boundary conditions are used. The 
exact solution is u{t,x) — r _rt *sin(.r — ct). We compute the solution up to T = 2, 
and use the LDG method with C defined by 


C = 




(2.9) 


We notice that, for this choice of fluxes, the approximation to the convective term 
cu x is the standard upwinding, and that the approximation to the diffusion term 
a d* u is the standard three point central difference, for the P° case. On the other 
hand, if one uses a central flux corresponding to e 12 = — c 2 i =0, one gets a spread- 
out five point central difference approximation to the diffusion term ad* u. 

The LDG methods based on with k = 1.2,3, 4 are tested. Elements with 
equal size are used. Time discretization is by the third-order accurate TVD Runge- 
Kutta method [38], with a. sufficiently small time step so that error in time is 
negligible comparing with spatial errors. We list the L 0 0 errors and numerical 
orders of accuracy, for as well as for its derivatives suitably scaled Ax m d™ Uh 
for 1 < /// < A*, at the center of of each element. This gives the complete description 
of the error for u over the whoh 1 domain, as in each element is a polynomial 
of degree A;. We also list the L ^ errors and numerical orders of accuracy for qh at 
the element center. 

In all the convection-diffusion runs with a > 0, accuracy of at least (A- + l)-th 
order is obtained, for both u /, and qh , when P k elements are used. See Tables 1 to 
3. The P 4 case for the purely convection equation a = 0 seems to be not in the as- 
ymptotic regime yet with N = 40 elements (further refinement with N = SO suffers 
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from round-off effects due to our choice of noil-orthogonal basis functions), Table 
4. However, the absolute values of the errors are comparable with the convection 
dominated case in Table 3. 


Table 1. The heat equation a = 1, c = 0. errors and numerical order of 

accuracy, measured at the center of each element, for Ax ni d" 1 u/, for 0 < m < k, 
and for <y/, . 


k 

variable 

N = 10 
error 

N = 20 

error order 

N = 40 

error order 


u 

4.55E-4 

5.79E-5 

2.97 

7.27E-6 

2.99 

1 

Ax d x u 

9.01E-3 

2.22E-3 

2.02 

5.5GE-4 

2.00 



4.17E-5 

2.48E-G 

4.07 

1.53E-7 

4.02 


u 

1.43E-4 

1.76E-5 

3.02 

2.19E-6 

3.01 

2 

Ax d T u 

7.87E-4 

1.03E-4 

2.93 

1.31E-5 

2.98 


(Ax) 2 d 2 u 

1.64E-3 

2.09E-4 

2.98 

2.62E-5 

2.99 


<i 

1.42E-4 

1.76E-5 

3.01 

2.19E-G 

3.01 


u 

1.54E-5 

9.6GE-7 

4.00 

G.11E-8 

3.98 


Ax d T u 

3.77E-5 

2.36E-6 

3.99 

1.47E-7 

4.00 

3 

(A.r) 2 d 2 u 

1.90E-4 

1.17E-5 

4.02 

7.34E-7 

3.99 


( A.r) ! dj u 

2.51E-4 

1.5GE-5 

4.00 

9.80E-7 

4.00 


Q 

1.48E-5 

9.66E-7 

3.93 

G.11E-8 

3.98 


a 

2.02E-7 

5.51E-9 

5.20 

1.63E-10 

5.07 


Ax dj -u 

1.65E-6 

5.14E-8 

5.00 

1.G1E-9 

5.00 

4 

(A.r) 2 d 2 u 

6.34E-G 

2.04E-7 

4.9G 

G.40E-9 

4.99 


( A.r) 3 O'l u 

2.92E-5 

9.47E-7 

4.95 

2.99E-8 

4.99 


(Ax) 4 d^u 

3.03E-5 

9.55E-7 

4.98 

2.99E-8 

5.00 


<1 

2.10E-7 

5.51E-9 

5.25 

1.G3E-10 

5.07 
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Table 2. The convection diffusion equation a = 1, c = 1. L lOC , errors and numerical 
order of accuracy, measured at the center of each element, for Aj* m d™ iih for 0 < 
m < A*, and for 


fc 

variable 

N = 10 
error 

N = 20 

error order 

N : 40 

error order 


u 

G.47E-4 

1.25E-4 

2.37 

1.59E-5 

2.97 

1 

Ax d r u 

9.G1E-3 

2.24E-3 

2.10 

5.5GE-4 

2.01 


<i 

2.9GE-3 

1.20E-4 

4.G3 

1.47E-5 

3.02 


a 

1.42E-4 

1.7GE-5 

3.02 

2.18E-6 

3.01 

2 

Ax dj it 

7.93E-4 

1.04E-4 

2.93 

1.31E-5 

2.99 


(A.r) 2 dj. u 

1.G1E-3 

2.09E-4 

2.94 

2.62E-5 

3.00 


< 1 

1.26E-4 

1.03E-5 

2.94 

2.12E-G 

2.95 


u 

1.53E-5 

9.75E-7 

3.98 

G.12E-8 

3.99 


Ax dj- u 

3.S4E-5 

2.34E-G 

4.04 

1.47E-7 

3.99 

3 

{Ax) 2 dju 

1.S9E-4 

1.18E-5 

4.00 

7.3GE-7 

4.00 


(A.; ) 3 d*u 

2.52E-4 

1.56E-5 

4.01 

9.81E-7 

3.99 


(1 

1.57E-5 

9.93E-7 

3.9S 

G.17E-8 

4.01 


ii 

2.04E-7 

5.50E-9 

5.22 

1.C4E-10 

5.07 


Ax d x u 

1.G8E-G 

5.19E-8 

5.01 

1.G1E-9 

5.01 

4 

( Ax ) 2 9 2 u 

G.3GE-6 

2.05E-7 

4.96 

G.42E-8 

5.00 


(Axf dj. u 

2.99E-5 

9.57E-7 

4.97 

2.99E-8 

5.00 


{Ax) 4 dju 

2.94E-5 

9.55E-7 

4.95 

3.00E-8 

4.99 


<1 

1.9CE-7 

5.35E-9 

5.19 

1.01E-10 

5.0G 
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Table 3. The convection dominated convection diffusion equation a = 0.01, c = 1. 
L 0 o errors and numerical order of accuracy, measured at the center of each element, 
for Ax m d™ ii h for 0 < tn < k, and for qh . 


k 

variable 

N = 10 
error 

N = 20 

error order 

N = 40 

error order 


u 

7.14E-3 

9.30E-4 

2.94 

1.17E-4 

2.98 

1 

A x djc U 

6.04E-2 

1.58E-2 

1.93 

4.02E-3 

1.98 


< i 

8.68E-4 

1.09E-4 

3.00 

1.31E-5 

3.05 


u 

9.59E-4 

1.25E-4 

2.94 

1.58E-5 

2.99 

2 

Ax d x u 

5.88E-3 

7.55E-4 

2.96 

9.47E-5 

3.00 


(Aj-) 2 d\u 

1.20E-2 

1.50E-3 

3.00 

1.90E-4 

2.98 


(l 

8.99E-5 

1.11E-5 

3.01 

1.10E-6 

3.34 


U 

1.11E-4 

7.07E-6 

3.97 

4.43E-7 

4.00 


Ax d x u 

2.52E-4 

1.71E-5 

3.88 

1.07E-6 

4.00 

3 

(Ax) 2 d^u 

1.37E-3 

8.54E-5 

4.00 

5.33E-6 

4.00 


(Axf dlu 

1.75E-3 

1.13E-4 

3.95 

7.11E-6 

3.99 


( l 

1.18E-5 

7.28E-7 

4.02 

4.75E-8 

3.94 


V 

1.85E-6 

4.02E-8 

5.53 

1.19E-9 

5.08 


Ax d r u 

1.29E-5 

3.76E-7 

5.10 

1.16E-S 

5.01 

4 

(Ax) 2 d‘ 2 u 

5.19E-5 

1.48E-6 

5.13 

4.65E-8 

4.99 


(Ax) 3 d 3 .u 

2.21E-4 

6.93E-6 

4.99 

2.17E-7 

5.00 


(Ax) 4 dju 

2.25E-4 

6.89E-6 

5.03 

2.17E-7 

4.99 


9 

3.58E-7 

3.06E-9 

6.87 

5.05E-11 

5.92 
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Table 4. The convection equation a . = 0, c — 1 . L 0 c errors and numerical order of 
accuracy, measured at the center of each element, for Ax m d™ iih for 0 < m < k. 


k 

variable 

N = 10 

N = 20 

N = 40 



error 

error 

order 

error 

order 

1 

u 

7.24E-3 

9.4GE-4 

2.94 

1.20E-4 

2.98 


Ax 0 r U 

G.09E-2 

1.60E-2 

1.92 

4.09E-3 

1.97 


u 

9.9GE-4 

1.28E-4 

2.9G 

1.G1E-5 

2.99 

2 

Ax d x u 

G.00E-3 

7.71E-4 

2.9G 

9.G7E-5 

3.00 


(A,-)' dl;u 

1.23E-2 

1.54E-3 

3.00 

1.94E-4 

2.99 


u 

1.2GE-4 

7.50E-6 

4.07 

4.54E-7 

4.05 

3 

Ax d x u 

1.G3E-4 

2.00E-5 

3.03 

1.07E-C 

4.21 



1.52E-3 

9.03E-5 

4.07 

5.45E-G 

4.05 


( A.r) ! d*u 

1.35E-3 

1.24E-4 

3.45 

7.19E-G 

4.10 


U 

3.55E-G 

8.59E-8 

5.37 

3.28E-10 

8.03 


Ax d x ii 

1.89E-5 

1.27E-7 

7.22 

1.54E-S 

3.05 

4 

(A.v) 2 dju 

8.49E-5 

2.28E-G 

5.22 

2.33E-S 

6.G1 


( A.r) :$ dju 

2.36E-4 

5.77E-G 

5.36 

2.34E-7 

4.G2 


(Ax) 4 d 4 a 

2.80E-4 

8.93E-G 

4.97 

1.70E-7 

5.72 


Finally, to show that the order of accuracy could really degenerate to k for P k \ 
as was already observed in [C], we rerun the heat equation case a . = Lc = 0 with 
the central flux 


C 


0 0 
0 0 


This time we can see that the global order of accuracy in L is only k when P k is 
used with an odd value of k. 
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Table 5. The heat equation a = 1, c = 0. L 0 0 errors and numerical order of 
accuracy, measured at the center of each element, for &x m d™ Uk for 0 < m < k , 


and for qh , using the central flux. 


k 

variable 

N = 10 

error 

N = 20 

error order 

N = 40 

error order 


a 

3.59E-3 

8.92E-4 

2.01 

2.25E-4 

1.98 

1 

Ax d T u 

2.10E-2 

1.06E-2 

0.98 

5.31E-3 

1.00 


(i 

2.39E-3 

6.19E-4 

1.95 

1.56E-4 

1.99 


a 

6.91E-5 

4.12E-G 

4.07 

2.57E-7 

4.00 

2 

A x d x ii 

7.CGE-4 

1.03E-4 

2.90 

1.30E-5 

2.98 


( A.r) 2 <9 2 u 

2.98E-4 

1.68E-5 

4.15 

1.03E-6 

4.02 


<1 

G.52E-5 

4.11E-G 

3.99 

2.57E-7 

4.00 


u 

1.62E-5 

1.01E-6 

4.00 

6.41E-8 

3.98 


Ax d T u 

1.0GE-4 

1.32E-5 

3.01 

1.G4E-6 


3 

( A.r) 2 dju 

1.99E-4 

1.22E-5 

4.03 

7.70E-7 

3.99 


(A.r) 3 dj u 

G.81E-4 

8.68E-5 

2.97 

1.09E-5 

2.99 


<2 

1.54E-5 

1.01E-6 

3.93 

G.41E-8 

3.98 


u 

8.25E-8 

1.31E-9 

5.97 

2.11E-11 

5.96 


A x d r a 

1.G2E-6 

5.12E-8 

4.98 

1.G0E-9 

5.00 

4 

(A.r) 2 dfu 

1.G1E-6 

2.41E-8 

G.0G 

3.78E-10 

G.00 


(A xf&'u 

2.90E-5 

9.46E-7 

4.94 

2.99E-8 

4.99 


(A.r) 4 d*u 

5.23E-6 

7.59E-8 

6.11 

1.18E-9 

6.01 


<1 

7.85E-8 

1.31E-9 

5.90 

2. 11E- 11 

5.96 


c. Proof of the nonlinear stability. In this section, we prove the the nonlinear 
stability result of Proposition 2.1. To do that, we first show how to obtain the 
corresponding stability result for the exact solution and then mimic the argument 
to obtain Proposition 2.1. 

The continuous case as a model. We start by rewriting the equations (2.5a) 
and (2.5b), in compact form. If in equations (2.5a) and (2.5b) we replace r u (.r) and 
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v q {x) by v u (x,t) and v ei (x,t). respectively, add the resulting equations, sum on j 


from 1 to A\ and integrate in time from 0 to T, we obtain that 


£?(w,v) = 0, V smooth v(/), V/e(0,7’), 


(2.10a) 


where 


#(w, v) = / f dtu(xj)v u (xj)dxdt+ f f q(x.t)v q (xj)d.rdt 
Jo Jo Jo Jo 

n h( w(.r, t ))* dj \(x.t)dx dt , 


(2.10b) 


using the fact that h(w(.r,/))* d r w(xj) = <9, ( o( u ) —g (u)q) is a complete deriva- 


tive, we see that 


B{ w, w) = 


f >r(x, 
Jo 


T ) dr + 


q-(xj)dxdt~- / u Q (x)dx , 
- Jo 


( 2 . 11 ) 


and that £>(w. w) = 0, by (2,10a), we immediately obtain the following Instability 


result: 


i f u 2 (x , t) dx + rr <f ( x . 

- Jo Jo Jo 


/)d.r J/ = - / w 0 (j-)dr. 


This is the argument we have to mimic in order to prove Proposition 2.1. 

The discrete case. Thus, we start by finding a compact form of equations 
(2.6a) and (2.6b). If we replace i'hdu(- r ) and by Vf^ u (x,t) and in 

the equations (2.6a) and (2.6b), add them up, sum on j from 1 to N and integrate 


in time from 0 to T, we obtain 


B h (w h . v*) = 0, Vv*«) € U* x 1/, V/ € (0, T). 


(2.12a) 


Bh('Wh'V h )=f f diUh(xJ)v h ,u{rJ)dx dt + f f <lh{x.t)vh t q(x,t) dx dt 
Jo Jo Jo Jo 

T 

- f V h(w*)J- +1 / 2 (f)[VA(/)]>+l / 2^ 

- '° !<><* 

- f V [ h(wh(xj))' dr Vh{x,t)dxdi. (2.1. 

k;<a- b 


(2.12b) 
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Next, we obtain an expression for <t?/,(w/,, W/,). It is contained in the following 
result. 

Lemma 2.3. We have 

B h {w h ,w h ) = 1 f u 2 h (x,T)dx f f q 2 h (xj)dxdt + 0r,c(K]) - 1 f u 2 h {x,0)dx, 

- Jo Jo Jo Jo 

where ©r,c([ w fe]) IS defined in Proposition 2.1. 

Next, since £J/,(w/, . w/, ) = 0, by (2.12a), we get the inequality 

\ [ u 2 h {x,T)dx+ f f ql(x,t)dxdt + 0 T .c([wft]) = ]- j u 2 h {x,0)dx 
- Jo Jo Jo ^ Jo 

from which Proposition 2.1 easily follows, since 

\ J u 2 h (.r,Q)dx < 1 J ul(x)dx , 

by (2.5c). It remains to prove Lemma 2.3. 

Proof Lemma 2.3. After setting V/, = W/, in (2.12b), we get 

B(w h ,w h ) =1 / u 2 h {x,T)dx+ f j ql(x,t)dxdt + f 9di,»{t) dt - 1 f ii 2 h (x,0)dx 
~ Jo Jo Jo Jo Jo 

where 

jh(w*)' +1/ . 2 (t) [w/,(t)] J+1/2 + f h(w fc (jr,*)) # 3r w A (ar,<)dar j. 

T 

It only remains to show that f 0 ®diss(t)dt = Or,c([ w /i])- 
To do that, we proceed as follows. Since 


h( ( X , f ) ) clj. W* ( iT , t ) ( f fll h ) ft ( ll ff ) Qh ) 'Uff ll f ) ) Of (Jh 

/«/ 1 

f(s)ds - g(u h )q h ) 

= d x (#«/,) -g(uh)qh) 

= d r H{w h (x,t)), 
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we get 


(-W*) = X | 

\ [77(wfe(/))] >+1/2 -h(w A )5 +1/2 (f)[wA(#)] i+ , /2 

1 <j<N 

l 

- E 1 

j [H(w h {t))} - h(wft)'(#)[w A (*)]| 

i<j<N 

L U+l/2 


Since, by the definition of H , 


[tf(w, ,(/))] = [ <i>(u h (t))} - [g(u h (t))q h {t)] 

= [c>(«a(*)) ] - - [</&(0 ]//(«/>(*))< 

and since (h u ,h q ) f = h, we get 

= X {[^("*(0)] - [^(»/i(0)]^/,(o - [«*(o] | 

1 <j<N ^ ' J+l/2 

+ X! {” [ 9 fc( # )]flf(“fc)( 0 -[</*(#)] ^? 

1<><N ^ 

77ns is the crucial step to obtain the L 2 -stability of the LDG methods , since the 
above expression gives us key information about the form that the flux h should 
have in order to make 0d 4 \«j,(t) a nonnegative quantity and hence enforce the In- 
stability of the LDG methods. Thus, by taking h as in (2.7a). we get 

@diss{t) = X f [ w *(^)] f( C [w*(#)] 
i<i<A ^ 

and the result follows. This completes the proof. □ 


} 


j+l/2 


} 


. 7 + 1/2 


This completes the proof of Proposition 2.1. 

c. The error estimate in the linear case. In this section, we prove the error 
estimate of Theorem 2.2 which holds for the linear case /'( ■ ) = c and a( • ) = a. To do 
that, we first show how to estimate the error between the solutions w„ = ( a u , (p , ) / . 
v — 1. 2, of 


d, u„ + d T (/(«„) - \Ja{u u )q v ) = 0 in (0.T) x (0. 1). 
q v - dj g{u„) = 0 in (0. T) x (0,1), 


u v {i = 0 ) = u 0iy . 


on (0. 1). 
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Then, we mimic the argument in order to prove Theorem 2.2. 

The continuous case as a model. By the definition of the form £(•,•), (2.10b), 
we have, for v = 1,2, 


£(w„,v) = 0, V smooth v(t), V# G (0,T). 

Since in this case, the form /?(•,•) is bilinear, from the above equation we obtain 
the so-called error equation : 


B(e,\) — 0. V smooth v(t), VtG(0,T), 


where e = Wj - w 2 . Now, from (2.11), we get that 

#( e 7 e )=^ [ e*(x,T)dz + f f e 2 q (x,t)dx dt [ c 2 u (x,0)dx, 
“ Jo Jo Jo J Jo 


and since e u (x,0) — «o,i( a ’) ~ i<o, 2 ( ;I ') and B(e,e) = 0, by the error equation , we 
immediately obtain the error estimate we sought: 



c u d x + 



e q {x,t) dx dt 


=-f 


( «o,i(*) - uoX x ) ) 2 dx - 


To prove Theorem 2.2, we only need to obtain a discrete version of this argument. 
The discrete case. Since. 


Bh{ w/,.v /f )=o, Vv A (t) g 14 x 14 , v#g(0,T), 

B h ( w,v/,) = o, Vv*(<) g 14 x Vh , v# g (o,T), 


by (2.12a) and by equations (2.5a) and (2.5b), respectively, we immediately obtain 
our error equation : 


B h {e,v h ) = 0, Vvft(f) G V* x 14, V/ G (0,T), 

where e = w — w . Now, according to the continuous case argument, we should con- 
sider next the quantity B h (e , e); however, since e is not in the finite element space, it 
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is more convenient to consider #/, (P/,(e), P/,(e)), where P*(e(f )) = ( P /,(<„(/))* P/,(e q [t)) ) 
is the so-called L 2 -projection of e(f) into the finite element space V£ x \ £ . The 
L 2 -projection of the function p into Vh , P h(p), is defined as the only element of the 
finite element space \) t such that 

f ( P k (p)(r) - p{ x ) ) v h (x)dx = 0, V v h e V h . (2.13) 

Jo 

Note that, in fact Uh(t — 0) = Fh(uo), by (2.Gc). 

Thus, by Lemma 2.3. we have 

1 T 1 

B h ( F h (e).F h (e))=l f \P h (<u(T))(x)\ 2 d,r+ f f \ F h (r q (t))(x) | 2 dx dt 
- Jo Jo Jo 

i r ] 

+ 0T,c([Pft(e)])-- / \F h (t u m(x)\ 2 dx. 

1 Jo 

and since 

Pft(e„(0)) = P*(«o - tu(0)) = Pft(tto) - Uh( 0) = 0, 
by (2.6c) and (2.13), and 

B h ( P*(e).P*(e))= 5*(Pfc(e)-e,P ft (e))= B h ( P fc (w) - w,P*(e)), 


by the error equation , we get 





(T))(.r) | 2 dx + 



+ 0T,c([P»(e)]) 


(t))(x)\ 2 dxdt 
- B h ( Pft(w) - w,Pft(e)). 


(2.14) 


Note that since in our continuous model, the right-hand side is zero, we expect the 
term £?(P/,(w) — w,P*(e)) to be small. 

Estimating the right-hand side. To show that this is so, we must suitably 
treat the term £?(f/,(w) — w, P^(e)). 
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Lemma 2.4. For p = P fc (w) - w, we have 

T 1 

5 A (p,P*(e)) < ^©r,c(p ) + \J J Q \F h (e q (t))(x)\ 2 dxdt 


+ (Ax) 


2k 


f 


C\(t) dt + ( Ax 




1/2 

r h (e u (t))(x)\ 2 dx } dt, 


wh 


ere 


CAO = 2 4 {( iia±£iiL Ar]ci2 l 2 dl )lu(t) ,2 + i +a(Aj .)2(i-*) N<) |2 + ] J, 

c 2 (t) = 2c fc | V^|ci2 |«(<)|*+2 +a(Ax) 2(i ‘ -fc) |«(<)|j [ +2 

where the constants c £ and dfc depend solely on k, and k = k except when the grids 
are uniform and k is even , in which case k = k + 1. 


Note how c n appears in the denominator of C\(t). However, C x (t) remains 
bounded as c\\ goes to zero since the convective numerical flux is an E-flux. 

To prove this result, we will need the following auxiliary lemmas. We denote by 
I u I//U-+ 1 )( j) the integral over J of the square of the (k + l)-the derivative of u. 

Lemma 2.5. For p == P/,(w) - w, we have 


I P u j+l/2 I ^ a-( A;r )^ + 1/2 | «. \ Hi k+i HJ . +1/s) , 

I [P« ]>+ 1/2 I < c k ( Ax )* +1 /- | u |//(fc + U( J i + 1/2 ), 

I Pij+ 1/2 I - c * V«( A;r )* +1/2 | « |//(*+ 2 ) ( j i+1/2) , 

I [P? b+l/2 I < \A* ( Ax ) k+i/ 2 | M | W iM-2)( J. + ]/2 ), 


where J J+ ]/ 2 — Ij U /j+i, t/ie constant c k depends solely on k, and k = A* except 
when the grids are uniform, and k is even, in which case k = k + 1. 

Proof. The two last inequalities follow from the first two and from the fact that 
q = s/o dj U. The two first inequalities with k = k follow from the definitions of 
and [ p u ] and from the following estimate: 


I Pft(«)(-* I i.j/ 2 ) ~~ U }+\/2 


< ^. k (Ax) k ^ 2 


u 


"'* +1 »(./,- +1/2 >* 
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where the constant c* depends solely on k. This inequality follows from the fact 
that P/i(* 0 ( ;r j+i /2 ) — u j+ 1/2 — 0 when t/ is a polynomial of degree k and from a 
simple application of the Bramble- Hilbert lemma; see Ciarlet [11]. 

To prove the inequalities in the case in which k = k + 1, we only need to show 
that if u is a polynomial of degree k + 1 for k even, then pff ~ 0. It is clear that it 
is enough to show this equality for the particular choice 

u(x) = ((.r-. r, + 1/2 )/(A.r/2)) A + I . 


To prove this, we recall that if P( denotes the Legendre polynomials of order (: 
(i) jl 1 Pf(- S ) P-m(*)d* = (ii ) iM±l) = (±1)', and (iii) P f (. s) is a linear 

combination of odd (even) powers of ,s for odd (even) values of ( . Since we are 
assuming that the grid is uniform, Axj = A.r y _|_j = Ax, we can write, by (i), 

P»(«)<*) = V ?l±i{ f PMuixi + 

o<f<k w w-i 

for x £ / p Hence, for our particular choice of t/. we have 
i of _i_ i f 1 

Kj+i, 2 = 2 E ^~ 2 ~j , &(*) a* - i ) k+1 p f (d + (.s + i)* +, p,(-i)}d* 

j 1 ) f Pr(s)* i {(-l) k+ '- i P t (l)+ Pr(-l)}ds 

f Pri*)* 1 {(-l) k+i -‘ +(-l) f }d« 

by (ii). When the factor { ( — 1)*+ 1-1 + ( — l) f } is different from zero, | k + 1 — / + f \ 
is even and since k is also even, | / — ( | is odd. In this case, by (iii), 

j Pf( s ) s 1 d.s = 0, 

and so pZj + — 0. This completes the proof. □ 

We will also need the following result that follows from a simple scaling argument; 


0 <(<k 

1 ^ 2f + 1 fk+V 

” 9 Z^ 9 

w ()<fj<k 

_ 1 +1 /A- + 1 

” 9 Z^ 9 


see Ciarlet [11]. 
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Lemma 2.6. We have 


I [P*(p)]i+i/2 I < 4 (Aar) -1 / 2 || P fc (p) || L2( j), 


where Jj+1/2 = Ij U h - M and the constant dk depends solely on k. 
We are now ready to prove Lemma 2.4. 


Proof of Lemma 2.4 ■ To simplify the notation, let us set v*, = P^e. By the defini- 
tion of £?/,(•,*), we have 


Bh(p<Vh) 


ff 


d t p u {xj)v llu (xj) dx (It + 


f T f P„(X, 
Jo Jo 


t) Vh,q{x,t) dx dt 


f ^(p) 5 +i/ 2 ( # )[va(*)]j+i/ 2* 

do i ^ ^ 


1<><A 7 


/ E / 


h(p(;r,f))* dj Vh{x J) dx dt 


T 

= -/ J ! Np )' + 1 / 2 (<)[ v ft (*)]_, +1/2 rff , 


l<j<A r 


by the definition of the L 2 -projection (2.13). 

Now, recalling that p = {puiPqY an< l that V/, = (v u ,i^)*\ we have 


h(p)' [v*(<)] = (cp^-cn [p u ])[<’„] 

+ (-y/ap^-c 12 [p,]) [t>« ] 
+ {-y/apZ + c\2 [p u ]) [v, ] 

= $1 + #2 + $3* 


By Lemmas 2.5 and 2.6, 

I 0i I < (A.r)*-' +1/2 | u | w *+i (<7) (|c| + cn) | [v u ] |, 

I 02 | < Q 4 ( A.r)*' (a | « |//A+2( jj ( A.r) A * + x/« I c i2 1 1 u |//*+ 2 (j))) || v u ||^2( 

I 03 I ^ c k dk ( Ax)* ( | « |//fe+i(J) (A^) fc k + | ('ll | | U |//fc + i(y)) ) || V q \\lH J). 
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This is the crucial step for obtaining our error estimates. Note that the treatment of 
9\ is very different than the treatment of #2 and 9%. The reason for this difference is 
that the upper bound for 9 X can be controlled by the form 0yvc([v/,])- we recall that 
V/, = IP/de). This is not the case for the upper bound for 9? because 0r,c[v/,] = 0 
if c = 0 nor it is the case for the upper bound for 9 3 because 0r.c[ v //] does not 
involve the jumps [v q ]l 

Thus, after a suitable application of* Young's inequality and simple algebraic 
manipulations, we get 


h(p)' [v A (i)] < [ v„ ] 2 + 1| 


I IMJ) 


+ ^C i (i)(Ax) 2k + ^C 2 (t)( Ax)' 


IU 2 (.7p 


Since 

Bh(p<Vh) < 

and since J jJr j j 2 — /, U Ij+ 1, the result follows after simple applications of the 
Cauchy- Schwartz inequality. This completes the proof. □ 


T 

[ Y I h (P)j+i/>(^) [ v * (#) li+1/2 \dt, 

JO 1 \ r 


1 <j<N 


Conclusion. Combining the equation (2.14) with the estimate of Lemma 2.4, 
we easily obtain, after a simple application of GronwalLs lemma, 

{/' IP .(..miWI’dr + jQ' | P k(c q (t))(x) I 2 dxdt + 0r.c([Pft(e)]) } ' 

< (A/)" jf v / C^)^ + (^’) A ' f,Mt |P /l ( f „(0)(.r)| 2 f/.r| ' dt. 


Theorem 2.2 follows easily from this inequality, Lemma 2.G, and from the following 
simple approximation result: 


||p-P/,(p) ||/.2 (0 ,i) < </*( Ax) k+] |p|//i«.+D(o,i) 


where g * depends solely on k\ see Ciarlet [11]. 
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3. The LDG methods for the multi-dimensional case. In this section, we 
consider the LDG methods for the following convection-diffusion model problem 

d t u + ^ d Ti (fi{u)- a ij(u) d Xj u) = 0 in (0, T) x (0, l) rf , (3.1a) 

\<i<d 1 <j<d 

u ( t = 0) = uo , on (0, l) rf , (3.1b) 

with periodic boundary conditions. Essentially, the one-dimensional case and the 
multidimensional case can be studied in exactly the same way. However, there are 
two important differences that deserve explicit discussion. The first is the treatment 
of the matrix of entries a jj(ii), which is assumed to be symmetric, semipositive 
definite and the introduction of the variables qc , and the second is the treatment of 
arbitrary meshes. 

To define the LDG method, w r e first notice that, since the matrix ci t j(u) is as- 
sumed to be symmetric and semipositive definite, there exists a symmetric matrix 
bij(u) such that 

dij(u) = bif(u) bf j{u). (3.2) 

l<t<d 

Then we define the new scalar variables qc — be ] (u) d Tj u and rewrite the 

problem (3.1) as follows: 

dt u + d Xi (/,(«) - bi({u)qt) = 0 in (0,T) x (0. l) d , (3.3a) 

i <i<d l <t<d 

qc- Y, d Xj g tj (u) = 0, ( = 1 in (0, T) x (0, l) d , (3.3b) 

1 <j<d 

u(t = 0) = u o , on (0,1 ^ (3.3c) 

where gcj{u) = f u bcj{s)ds. The LDG method is now obtained by discretizing 
(3.3) by the Discontinuous Galerkin method. 

We follow what was clone in §2. So, we set w = (u, q) f = (u,q\,-- • , q^Y and, 
for each i — 1, • • * , c/, introduce the flux 

hi(w )=(/«(«)- L b iti u )qt, -gu(u),' • ■ , - 9 di{u ) )'. (3.4) 

1 <f <d 
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We consider triangulations of (0, l) d , Taj — { K }, made of non-overlapping poly- 
hedra. We require that for any two elements K and K\ K fl K* is either a face 
e of both A and A' with nonzero (d — 1)-Lebesgue measure | e |, or has Hausdorff 
dimension less than d — 1. We denote by Sax the set of all faces e of the border 
of A for all K £ T“at • The diameter of A is denoted by A.r/c and the maximum 
A./* a*, for A £ 7 a j is denoted by Ax. We require, for the sake of simplicity, that 
the triangulations Tax be regular, that is, there is a constant independent of A.r 
such that 

— < tr VA € r AJ . 

PK 

where pj^ denotes the diameter of the maximum ball included in A . 

We seek an approximation W/, — (u/,,q/,) / = (?/& , g/,i ■> • * * . ^/^/ ) / to w such that 
for each time t £ [0. T}. each of the components of w/, belong to the finite^ element 
space 

Vh = v h k = { v€ !'(((), 1 ) d ) : r\ h G P k '{K) V A' € Ta,}. (3.5) 


where' P*(A~) denotes the space of polynomials of total degree at most A*. In or- 
der to determine the approximate solution W/ M we proceed exactly as in the one- 
dimensional case. This time, however, the integrals are made on each element A of 
the triangulation 7 a j - We obtain the following weak formulation on each element 
A^ of the triangulation 7 a t* 


/ d f u h {xj) v h , a (x)dx - V* / /i rti (w /? (.rd))c), ( . v hM (x)dx 

Jl < X 

+ f /»u(wfc,iifl/ v -)(:r,/) v h ' U {x)dT{.r) = 0, V r/,, u G P k (I\), 
JciK 


q h f(.rj)v, hqi (x)dx- ^ (.r)r/.r 

+ l hqc( w h^nnh){xj) v h ^ f (x)dT(x) = 0, V G P k (I\). 
Jdl< 


(3.6a) 


(3.6b) 
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L 

where vidK denotes the outward unit normal to the element K at x (E dK. It 
remains to choose the numerical flux {h u ,h. qi , • ■ ■ ,hq d Y = h = h(w/, , naA')(^^)- 
As in the one- dimensional case, we require that the fluxes h be of the form 

h(wft , n#A-)(.r) = h(wk(x tn1h \t).w h (x exiK ,t);n dK ), 

where w^(x iniR ) is the limit at x taken from the interior of K and w*(a: cr<A ') the 
limit at x from the exterior of A , and consider fluxes that (i) are locally Lipschitz, 
conservative, that is, 

h(w h (x inl “ ),w k (x e * tK ); n aK ) + h(w h (x extK ),w h (x intK ); -n aK ) = 0, 


<//, (j-.O ) d. 


• T = / 
Jk 


uq(x) v h ,i(x)dx , V v h j € P fc ( A'), 


(3.6c) 


and consistent with the flux 

^ ^ hi W'r)/v\i, 
l<i<d 

(ii) allow for a local resolution of each component of q/, in terms of tih only , (iii) 
reduce to an E-flux when «(•) = 0, and that (iv) enforce the L 2 -stability of the 
method. 

Again, we write our numerical flux as the sum of a convective flux and a diffusive 
flux: 


h — h COTM? + h difff 

where the convective flux is given by 

h C o»«.(w _ , w + ; n) = (/(«”, u + ; n),0)', 

where f( u~ , u + ; n) is any locally Lipschitz E-flux which is conservative and consis- 
tent, with the nonlinearity 

53 /«(»)»«• 

!<«<(/ 
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and the diffusive flux h<n//(w , w + ; n) is given by 


(- 


E 

1 <ij<d 


!»«/(«)] 


qt n i‘ 


\<i<d 


X 9 id(u)n, )' - Cdiff [w], 

1<T<rf 


where 


Cdiff 


/ 0 c i2 ej3 

j -r 12 0 0 

-rn 0 0 


fi<t\ 

0 

0 


\-cid 0 0 ••• 0 / 


c\j = ri j( w . w+ ) is locally Lipschitz for j ~ 1, ■ * • . (U 


c\j = 0 when r/(-) = 0 for j = 1. • ■ • , d. 


We claim that this flux satisfies the properties (i) to (iv). 

To prove that properties (i) to (iii) are satisfied is now a simple* exercise. To see 
that the property (iv) is satisfied, we first rewrite the flux h in the following way: 


(- E 


[<7if(*0] — V" 7 — 7 

r l (K Hi. - > ffi 1 ( U ) V{, ■ 

i/ 

KiJ<d L J l<i<d 


, - ^2 n * y - c[ 


w 


I <i<d 


where 


c 


/ ('ll Cl-2 Ci3 ■■■ 

—c\> 0 0 • • ■ 

-r, ; , 0 0 ••• 


<‘\d \ 
0 
0 


\-c,rf 0 0 ••• 0 7 


cii = 


h( E 


Pi\U 


a I \ *■ — 4 f u 1 
i < i < rf 1 J 


n t — f(u ,f/ + ;n) 


where o,( u) = f u /,(.s) <7.s. Since /(-,-;n) is an E-flux, 

1/‘ M + 

cii = j— [7 / ( XI /'(*)”« _ /(■«“, « + ;n) ) ds > 0. 

l^J ./»- i<i<rf 


and so the matrix C is semipositive definite. The property (iv) follows from this 
fact and from the following multidimensional version of Proposition 2.1. 



30 


Proposition 3.1. (L 2 -stabilitv) We have, 


\ f u 2 h {x,T)dx+ f f \q h (xj)\ 2 dxdt + 6 T ,£{[wh]) <\ f ul{x)dx , 
- ^(0,l) d JO J{ 0,1) J - J( 0,l) d 


where 


T 

0r.d[w h })= I V f [w h (xJ)} t C[w h (xJ)]dT(x)dt. 

Jo 

We can also prove the following error estimate. We denote the integral over 
(0,l) rf of the sum of the squares of all the derivatives of order (k + 1) of it by 


[2 

T + 1 ‘ 


Theorem 3.2. (L 2 -error estimate) Let e be the approximation error w — w/, . Then 
we have , for arbitrary , regular grids , 

| f | e u (x, T) | 2 dx + f f | e q (x,t) | 2 dxdt + Qr,e([e])} < C(Aa-) fc , 

L J(0,l) d Jo J( 0A) d J 

where C = C'{h\ \ u | u | a -+-2 ) • In the purely hyperbolic case a i] = 0. the constant 
C is of order (Aa*) 1 / 2 . In the purely parabolic case c — 0 T the constant C is of order 
Ax for even values of k and of order 1 otherwise for Cartesian products of uniform, 
grids and for C identically zero provided that the local spaces Q k are used instead 
of the spaces P k , where Q k is the space of tensor products of one dimensional 
polynomials of degree k. 


The algebraic manipulations needed to prove this result are a straightforward 
extension to the multidimensional case of the manipulations of the proof of the 
corresponding one- dimensional result, Theorem 2.2. The approximation properties 
of the finite element spaces Vh that extend the results of Lemmas 2.5 and 2.6 are 
the following. Let e denote a face of the element K and let us denote by K e the 
element sharing the face e with A’, then 

|| P/j(a) ± - u \\l 2 (€) < ( tAx ) A + 1 / 2 | u |tf(*+i)(/ v 'uA^p 
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where Fh.(u)± is either the value of IP/, ( 1 / ) at e from the interior of K or from its 
exterior, and 


II [P/i(p)] \\l 2 (c) < dk{Ax) || Fh(p) ||l 2 (A'uA\,p 

where [P/ J ( J r>)] denotes the jump at c of P h{p)- Finally, we also use the following 
result: 

II p — ^ h{p) ||L 2 (0,l) d — ffk { A*' ) k + ] I P |//<fr + O( 0i i)rf, 

All these approximation results can he found in Ciarlet [11], for example. 

4 . Concluding remarks. I 11 this paper, we have considered the so-called LDG 
methods for convection-diffusion problems. For scalar problems in multidimensions, 
we have shown that they are L 2 -stable and that in the linear case, they are of order 
k if polynomials of order k are tised. We have also shown that this estimate is sharp 
and have displayed the strong dependence of the order of convergence of the LDG 
methods on the choice of the numerical fluxes. 

The LDG methods for multidimensional systems, like for example the compress- 
ible Navier- Stokes equations and the equations of the hydrodynamic model for 
semiconductor device simulation, can be easily defined by simply applying the pro- 
cedure described for the multidimensional scalar case to each component of u. In 
practice, especially for viscous terms which are not symmetric but still semipos- 
itive definite, such as for the compressible Navier-Stokes equations, we can use 
q — 3 j f/ u) as the auxiliary variables. Although with this choice, the L 2 - 

stability result will not bo available theoretically, this would not cause any problem 
in practical implementation; see Bassi and Rebay [5] and Bassi et al [6]. 

The main advantage of these methods is their extremely high parallelizabil- 
ity and their high-order accuracy which render them suitable for computations 
of convection-dominated flows. Indeed, although the LDG method have a large 
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amount of degrees of freedom per element, and hence more computations per ele- 
ment are necessary, its extremely local domain of dependency allows a very efficient 
parallelization that by far compensates for the extra amount of local computations. 

Acknowledgments. The authors want to thank X. Makridakis for bringing to 
theii attention the reference [3], and H. Atkins for helpful discussions in the com- 
putation. 
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